R Notebook

source("include.R")

❤ Loaded agrimonia dataset. Available as df_agri.
❤ Converted df_agri$Time to Date variable type, year-month-day.
❤ Created DATE_FORMAT variable for date comparisons.
❤ Loaded stations split dataset. Available as df_stat.

Maps plots

library(maps)

Avvertimento: il pacchetto ‘maps’ è stato creato con R versione 4.2.3

library(ggplot2)

Caricamento pacchetto: ‘ggplot2’

Il seguente oggetto è mascherato da ‘package:crayon’:

    %+%

library(sf)

Avvertimento: il pacchetto ‘sf’ è stato creato con R versione 4.2.3Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

Idea: plot the values for pm10 each day, like with “bubbles” of ggplot, so to get 2192 plots. And then combine all them in a single gif, which shows the evolution over time of those pm10 values.

https://askubuntu.com/questions/648244/how-do-i-create-an-animated-gif-from-still-images-preferably-with-the-command-l

So idea: plot bubbles with size according to the pm10 value. And maybe if we want use also a different color (like to separate high/medium/low values).

regioni_italiane <- st_read("italia/gadm40_ITA_1.shp") 

Reading layer `gadm40_ITA_1' from data source 
  `C:\Users\feder\Desktop\Uni magistrale\Bayesian statistics\progetto-bayesian\src\italia\gadm40_ITA_1.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 20 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 6.630879 ymin: 35.49292 xmax: 18.52069 ymax: 47.09096
Geodetic CRS:  WGS 84

sites <- data.frame(
    longitude = unique(df_agri$Longitude), 
    latitude = unique(df_agri$Latitude))

pad = 0.2*c(-1,1)
mappa_completa <- ggplot() +
  geom_sf(data = regioni_italiane, fill = "lightblue", color = "black", size = 0.5)+
      coord_sf(xlim = range(sites$longitude)+pad, 
             ylim = range(sites$latitude)+pad, expand = FALSE)+
  # geom_sf(data = world) +
  geom_point(data = sites, aes(x = longitude, y = latitude), size = 1, shape = 23, fill = "darkred")
 
print(mappa_completa)

sites <- data.frame(
    longitude = unique(df_agri$Longitude), 
    latitude = unique(df_agri$Latitude))

for (data in seq.Date(from=as.Date("2016-01-01",DATE_FORMAT),
               to=as.Date("2016-01-11",DATE_FORMAT),by="days")){
        
    dt = as.Date(data,origin="1970-01-01")
    df_date = df_agri[which(df_agri$Time == dt),]
    print(dt)
    
    pad = 0.2*c(-1,1)
    mappa_completa <- ggplot() +
      geom_sf(data = regioni_italiane, fill = "lightblue", color = "black", size = 0.5)+
          coord_sf(xlim = range(sites$longitude)+pad, 
                 ylim = range(sites$latitude)+pad, expand = FALSE)+
      # geom_sf(data = world) +
      geom_point(data = sites, aes(x = longitude, y = latitude), size = 1, shape = 23, fill = "darkred")+
        geom_point(data=df_date,aes(x=Longitude,y=Latitude),size=df_date$AQ_pm10/10,alpha=0.5)+
        ggtitle(dt)
    
    print(mappa_completa)
    ggsave(paste0("./maps_images/",dt,".jpeg"),plot=mappa_completa)
}

[1] "2016-01-01"
[1] "2016-01-02"
[1] "2016-01-03"
[1] "2016-01-04"
[1] "2016-01-05"
[1] "2016-01-06"
[1] "2016-01-07"
[1] "2016-01-08"
[1] "2016-01-09"
[1] "2016-01-10"
[1] "2016-01-11"

Trend plots

stations = unique(df_agri$IDStations)
cols = colora(6,97,show=1)

for (st in stations[1:4]){
    df_st = df_stat[[st]]
    data2016 = subset(df_st,Time>="2016-01-01" & Time<="2016-12-31")
    data2017 = subset(df_st,Time>="2017-01-01" & Time<="2017-12-31")
    data2018 = subset(df_st,Time>="2018-01-01" & Time<="2018-12-31")
    data2019 = subset(df_st,Time>="2019-01-01" & Time<="2019-12-31")
    data2020 = subset(df_st,Time>="2020-01-01" & Time<="2020-12-31")
    data2021 = subset(df_st,Time>="2021-01-01" & Time<="2021-12-31")
    
    
    plot(as.Date(format(data2016$Time,"%d %m"),"%d %m"),data2016$AQ_pm10 ,type="l",col=cols[1],
         xlab="time",ylab="PM_10 values",main=paste0("Station ",st,", all years"),
         ylim=range(na.omit(df_st$AQ_pm10)))
    lines(as.Date(format(data2017$Time,"%d %m"),"%d %m"),data2017$AQ_pm10,type="l",col=cols[2])
    lines(as.Date(format(data2018$Time,"%d %m"),"%d %m"),data2018$AQ_pm10,type="l",col=cols[3])
    lines(as.Date(format(data2019$Time,"%d %m"),"%d %m"),data2019$AQ_pm10,type="l",col=cols[4])
    lines(as.Date(format(data2020$Time,"%d %m"),"%d %m"),data2020$AQ_pm10,type="l",col=cols[5])
    lines(as.Date(format(data2021$Time,"%d %m"),"%d %m"),data2021$AQ_pm10,type="l",col=cols[6])
    legend('topright', levels(as.factor(2016:2021)), fill=cols, bty='n', cex = 0.6)

}

cols = colora(10,show=0)
for (year in 2016:2021){
    for (i in 1:length(stations)){
        st = stations[i]
        df_st = df_stat[[st]]
        df_st = df_st[df_st$Time>=as.Date(paste0(year,"-01-01"),"%Y-%m-%d") &
                      df_st$Time<=as.Date(paste0(year,"-12-31"),"%Y-%m-%d"),]
        if(st=="1264"){
            plot(df_st$Time,df_st$AQ_pm10,type="l",col=cols[i%%10],
                 ylim=range(na.omit(df_agri$AQ_pm10)),main=paste0(year," all the stations"),
                 xlab="months",ylab="PM_10 values")
        } else{
            lines(df_st$Time,df_st$AQ_pm10,type="l",col=cols[i%%10])
        }
    }
}